% 双线性变换法
close all;
clear;
clc;
%给定数字滤波器指标
omega_p = 0.2*pi;            %通带截止频率
omega_st = 0.3*pi;            %阻带截止频率
R_p = 1       ;   %通带允许的最大衰减
R_st = 15      ;   %阻带允许的最小衰减
f_s = 1000    ;   %采样频率
T_s = 1/f_s ;   %采样间隔

%预畸变处理,将归一化数字角频率参数变换成模拟低通滤波器的角频率参数
C = 2/T_s ;
Omega_p = C * tan(omega_p/2 );
Omega_st = C * tan(omega_st/2 ); 
%对模拟低通滤波器的角频率参数做归一化处理
lamda_p  = Omega_p  / Omega_p;
lamda_st = Omega_st / Omega_p;
%设计归一化模拟滤波器,得到归一化模拟低通滤波器的转移函数
[N,Wn] = buttord( Omega_p , Omega_st , R_p , R_st , 's' ); %选择模拟巴特沃斯低通滤波器的最小阶数
[z,p,k] = buttap(N); %创建巴特沃斯模拟低通滤波器
[Bp,Ap] = zp2tf(z,p,k); %由零点、极点、增益确定传输函数的分子与分母系数
%利用归一化模拟低通滤波器的转移函数确定模拟低通滤波器的转移函数
[b,a] = lp2lp(Bp,Ap,Wn);
%利用模拟低通滤波器的转移函数确定IIR数字滤波器的转移函数 
[bz,az] = bilinear(b,a,f_s);

figure(1);
freqz(bz,az);
title('低通滤波器幅度谱和相位谱特性');

%滤波效果测试 
N = 1000;
t = [0:N-1]*T_s ;
f1 = 50;
f2 = 250;
x1 = sin(2*pi*f1*t);
x2 = sin(2*pi*f2*t);
x  = x1 + x2;

fft_x = fft(x);
X_mag = fftshift(abs(fft_x))/N ;
X_ang = fftshift(angle(fft_x));
delta_f = f_s/N;
f = ( -N/2:N/2-1)*delta_f;

figure(2);
plot(t,x);
title('原信号时域波形');
xlabel('t(s)');

%滤波
lp_x = filter(bz,az,x);
lp_fft_x = fft(lp_x);
lp_X_mag = fftshift(abs(lp_fft_x))/N ;
lp_X_ang = fftshift(angle(lp_fft_x));
figure(3);
plot(t,lp_x);
title('滤波后信号时域波形');
xlabel('t(s)');
